﻿#define _CRT_SECURE_NO_WARNINGS
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <windows.h>
#include "BasicMath.h"
#include "Shlwapi.h"
#pragma comment(lib,"Shlwapi.lib") 


/*实现MatBase中方法*/
MatBase::MatBase() 
{
	OnInit = 0;
	rows = 0;
	cols = 0;
	data = NULL;
	//related = false;
}
MatBase::MatBase(unsigned int _rows, unsigned int _cols, int _OnInit)
{
	OnInit = _OnInit;
	rows = _rows;
	cols = _cols;
	if (_OnInit)
		data = (real*)malloc(rows * cols * sizeof(real));
	//related = false;
}
MatBase::MatBase(const MatBase& mat)
{
	OnInit = 1;
	rows = mat.rows;
	cols = mat.cols;
	if (mat.OnInit == 1) //复制永久对象
	{
		data = (real*)malloc(mat.rows * mat.cols * sizeof(real));
		memcpy(data, mat.data, rows * cols * sizeof(real));
	}
	else if (mat.OnInit == 2)  //传递地址，销毁临时对象
	{
		data = mat.data;
	}
	else
		_ERROR("对象赋值", "赋值右值未初始化");
}
void MatBase::Init(unsigned int _rows, unsigned int _cols, int _OnInit)
{
	OnInit = _OnInit;
	rows = _rows;
	cols = _cols;
	if (_OnInit)
		data = (real*)malloc(rows * cols * sizeof(real));
	//related = false;
}
void MatBase::Delete() 
{
	if (OnInit)
	{
		OnInit = 0;
		rows = 0;
		cols = 0;
		if (data)
			free(data);
		data = NULL;
	}
}
unsigned int MatBase::Rows() const
{
	return rows;
}
unsigned int MatBase::Cols() const
{
	return cols;
}
real* MatBase::PointElem(unsigned int x, unsigned int y) const
{
	return POINTELEM(data, cols, x, y);
}
real MatBase::GetElem(unsigned int x, unsigned int y) const
{
	if (!CheckLimit1(x, 1, rows) || !CheckLimit1(y, 1, cols))
	{
		_ERROR("取矩阵元素", "矩阵下标超出范围");
		return 0;
	}
	return *POINTELEM(data, cols, x, y);
}
bool MatBase::SetElem(unsigned int x, unsigned int y, const real val)
{
	if (!CheckLimit1(x, 1, rows) || !CheckLimit1(y, 1, cols)) 
	{
		_ERROR("矩阵元素赋值", "矩阵下标超出范围");
		return false;
	}
	*POINTELEM(data, cols, x, y) = val;
	return true;
}
void MatBase::SetAll(const real val)
{
	for (int i = 0;i < rows * cols;i++)
		data[i] = val;
}
void MatBase::SetZeros()
{
	memset(data, 0, rows * cols * sizeof(real));
}
/*
void MatBase::Copy(const MatBase& mat) 
{
	if (rows != mat.rows || cols != mat.cols)
		_ERROR("矩阵赋值", "行列数不匹配");
	memcpy(data, mat.data, rows * cols * sizeof(real));
}
*/
void MatBase::Print() const
{
	if (!OnInit) 
	{
		_ERROR("输出矩阵","未初始化");
		return;
	}
	for (int i = 1;i <= rows;i++) 
	{
		for (int j = 1;j <= cols;j++)
			SHOWREAL(*POINTELEM(data, cols, i, j));
		ENDL;
	}
}
MatBase MatBase::GetTrans() const
{
	MatBase ans(cols, rows);
	for (int x = 1;x <= rows;x++)
		for (int y = 1;y <= cols;y++)
			*POINTELEM(ans.data, ans.cols, y, x) = *POINTELEM(data, cols, x, y);
	ans.OnInit = 2;
	return ans;
}
bool MatBase::SetTrans() 
{
	if (!OnInit) 
	{
		_ERROR("矩阵转置","未初始化");
		return false;
	}
	real* temp = (real*)malloc(rows * cols * sizeof(real));
	for (int x = 1;x <= rows;x++)
		for (int y = 1;y <= cols;y++)
			*POINTELEM(temp, rows, y, x) = *POINTELEM(data, cols, x, y);
	free(data);
	data = (real*)malloc(rows * cols * sizeof(real));
	memcpy(data, temp, rows * cols * sizeof(real));
	free(temp);
	unsigned int t = rows;
	rows = cols;
	cols = t;
	return true;
}
MatBase MatBase::GetExpand(int _rows, int _cols, int _x, int _y, real fill) const
{
	if (_x + rows - 1 > _rows || _y + cols - 1 > _cols)
	{
		_ERROR("矩阵扩展","行列数过小，矩阵溢出");
		return MatBase();
	}
	MatBase ans(_rows, _cols);
	if (fill)
		ans.SetAll(fill);
	else
		ans.SetZeros();
	for (int i = 1;i <= rows;i++)
		memcpy(POINTELEM(ans.data, ans.cols, i + _x - 1, _y), POINTROW(data, cols, i), cols * sizeof(real));
	ans.OnInit = 2;
	return ans;
}
bool MatBase::SetExpand(int _rows, int _cols, int _x, int _y, real fill)
{
	if (_x + rows - 1 > _rows || _y + cols - 1 > _cols)
	{
		_ERROR("矩阵扩展", "行列数过小，矩阵溢出");
		return false;
	}
	real* tempData = (real*)malloc(rows * cols * sizeof(real));
	unsigned int tempRows = rows;
	unsigned int tempCols = cols;
	memcpy(tempData, data, rows * cols * sizeof(real));
	free(data);
	data = (real*)malloc(_rows * _cols * sizeof(real));
	rows = _rows;
	cols = _cols;
	if (fill)
		SetAll(fill);
	else
		SetZeros();
	for (int i = 1;i <= tempRows;i++)
		memcpy(POINTELEM(data, cols, i + _x - 1, _y), POINTROW(tempData, tempCols, i), tempCols * sizeof(real));
	free(tempData);
	return true;
}
bool MatBase::SaveBIN(const char* szPath)
{
	if (!OnInit)
	{
		_ERROR("保存二进制文件", "矩阵未初始化");
		return false;
	}
	FILE* file;
	file = fopen(szPath, "wb");
	if (file == NULL) 
	{
		_ERROR("保存二进制文件","无法打开文件");
		return false;
	}
	fwrite(&rows, sizeof(unsigned int), 1, file);
	fwrite(&cols, sizeof(unsigned int), 1, file);
	fwrite(data, rows * cols * sizeof(real), 1, file);
	fclose(file);
	return true;
}
bool MatBase::LoadBIN(const char* szPath)
{
	FILE* file;
	file = fopen(szPath, "rb");
	if (file == NULL)
	{
		_ERROR("读取二进制文件", "无法打开文件");
		return false;
	}
	if (OnInit)
		free(data);
	else
		OnInit = 1;
	fread(&rows, sizeof(unsigned int), 1, file);
	fread(&cols, sizeof(unsigned int), 1, file);
	data = (real*)malloc(rows * cols * sizeof(real));
	fread(data, rows * cols * sizeof(real), 1, file);
	fclose(file);
	return true;
}
bool MatBase::SaveTXT(const char* szPath) 
{
	if (!OnInit)
	{
		_ERROR("保存TXT文件", "矩阵未初始化");
		return false;
	}
	FILE* file;
	file = fopen(szPath, "w");
	if (file == NULL)
	{
		_ERROR("保存TXT文件", "无法打开文件");
		return false;
	}
	for (int i = 1;i <= rows;i++) 
	{
		for (int j = 1;j <= cols;j++)
			fprintf(file, "%.12f ", *POINTELEM(data, cols, i, j));
		fprintf(file, "\n");
	}
	fclose(file);
	return true;
}

bool MatBase::LoadTXT(const char* szPath)
{
	FILE* file;
	file = fopen(szPath, "r");
	if (file == NULL)
	{
		_ERROR("读取TXT文件", "无法打开文件");
		return false;
	}
	if (OnInit)
		free(data);
	else
		OnInit = 1;
	rows = 0;
	int c0 = -1, c1;
	char* buffer = (char*)malloc(MAXROW * sizeof(char));
	while (fgets(buffer, MAXROW * sizeof(char), file) != NULL)
	{
		rows++;
		c1 = 0;
		for (int i = 0;i < MAXROW * sizeof(char) - 1;i++)
		{
			if (buffer[i + 1] == ' ' || buffer[i + 1] == '\t')
			{
				if (buffer[i] != ' ' && buffer[i] != '\t')
					c1++;
			}
			else if (buffer[i + 1] == '\r' || buffer[i + 1] == '\n' || buffer[i + 1] == '\0') 
			{
				if (buffer[i] != ' ' && buffer[i] != '\t')
					c1++;
				break;
			}
		}
		if (c0 != -1 && c0 != c1) 
		{
			_ERROR("读取TXT文件", "每行数据列数不一致");
			OnInit = 0;
			return false;
		}
		c0 = c1;
	}
	free(buffer);
	cols = c0;
	data = (real*)malloc(rows * cols * sizeof(real));
	rewind(file);
	for (int i = 1;i <= rows;i++)
		for (int j = 1;j <= cols;j++)
			fscanf(file, "%lf", POINTELEM(data, cols, i, j));
	fclose(file);
	return true;
}
bool MatBase::SaveCSV(const char* szPath)
{
	if (!OnInit)
	{
		_ERROR("保存CSV文件", "矩阵未初始化");
		return false;
	}
	FILE* file;
	file = fopen(szPath, "w");
	if (file == NULL)
	{
		_ERROR("保存CSV文件", "无法打开文件");
		return false;
	}
	for (int i = 1;i <= rows;i++)
	{
		for (int j = 1;j <= cols;j++)
			if (j == cols)
				fprintf(file, "%.12f\n", *POINTELEM(data, cols, i, j));
			else
				fprintf(file, "%.12f,", *POINTELEM(data, cols, i, j));
	}
	fclose(file);
	return true;
}

bool MatBase::LoadCSV(const char* szPath)
{
	FILE* file;
	file = fopen(szPath, "r");
	if (file == NULL)
	{
		_ERROR("读取CSV文件", "无法打开文件");
		return false;
	}
	if (OnInit)
		free(data);
	else
		OnInit = 1;
	rows = 0;
	int c0 = -1, c1;
	char* buffer = (char*)malloc(MAXROW * sizeof(char));
	while (fgets(buffer, MAXROW * sizeof(char), file) != NULL)
	{
		rows++;
		c1 = 0;
		for (int i = 0;i < MAXROW * sizeof(char);i++)
		{
			if (buffer[i] == '\r' || buffer[i] == '\n' || buffer[i] == '\0')
				break;
			if (buffer[i] == ',')
				c1++;
		}
		if (c0 != -1 && c0 != c1)
		{
			_ERROR("读取CSV文件", "每行数据列数不一致");
			OnInit = 0;
			return false;
		}
		c0 = c1;
	}
	free(buffer);
	cols = c0 + 1;
	data = (real*)malloc(rows * cols * sizeof(real));
	rewind(file);
	for (int i = 1;i <= rows;i++)
		for (int j = 1;j <= cols;j++)
			if (j == cols)
				fscanf(file, "%lf", POINTELEM(data, cols, i, j));
			else
				fscanf(file, "%lf,", POINTELEM(data, cols, i, j));
	fclose(file);
	return true;
}

void MatBase::Plot(unsigned int type, unsigned int index, const char* tittle)
{
	//判断异常
	if (!OnInit)
	{
		_ERROR("绘制图像", "未初始化");
		return;
	}
	if (!PathFileExistsA(PATH_PLOT))
	{
		char szBuffer[64];
		sprintf(szBuffer, "找不到画图模块 %s", PATH_PLOT);
		_ERROR("绘制图像", szBuffer);
		return;
	}
	if ((type == ROW && index > rows) || (type == COL && index > cols))
	{
		_ERROR("绘制图像", "坐标超限");
		return;
	}
	//采样压缩
	unsigned int length0 = (type == ROW) ? cols : rows;
	unsigned int sumb = GetBitNum(length0);
	unsigned int ite = 1;
	if (sumb > 3)
		ite = pow(10, sumb - 3.0);
	unsigned int length = length0 / ite;
	//进程通信管道
	SECURITY_ATTRIBUTES saAttr;
	saAttr.nLength = sizeof(SECURITY_ATTRIBUTES);
	saAttr.bInheritHandle = TRUE;
	saAttr.lpSecurityDescriptor = NULL;
	HANDLE handle_read;
	HANDLE handle_write;
	bool ret = CreatePipe(&handle_read, &handle_write, &saAttr, 0);
	if (!ret)
	{
		_ERROR("绘制图像", "创建进程通信管道失败");
		return;
	}
	//启动绘图模块进程
	char cmdline[] = PATH_PLOT;
	PROCESS_INFORMATION piProcInfo;
	ZeroMemory(&piProcInfo, sizeof(PROCESS_INFORMATION));
	STARTUPINFO si;
	ZeroMemory(&si, sizeof(STARTUPINFO));
	si.cb = sizeof(STARTUPINFO);
	si.hStdInput = handle_read; //把管道的读句柄传给子进程
	si.dwFlags |= STARTF_USESTDHANDLES;
	ret = CreateProcess(NULL, cmdline, NULL, NULL, TRUE, CREATE_NEW_CONSOLE, NULL, NULL, &si, &piProcInfo);
	if (!ret) 
	{
		_ERROR("绘制图像", "启动绘图模块进程失败");
		return;
	}
	else
	{
		CloseHandle(piProcInfo.hProcess);
		CloseHandle(piProcInfo.hThread);
		CloseHandle(handle_read);
	}

	//写入文件头
	FigHead fh;
	memcpy(fh.szTitle, tittle, (strlen(tittle) + 1) * sizeof(char));
	fh.length = length;
	DWORD len;
	ret = WriteFile(handle_write, &fh, sizeof(FigHead), &len, NULL); //子进程读了后，父进程才可以继续写入管道
	if (!ret)
	{
		_ERROR("绘制图像", "写入文件头失败");
		return;
	}

	//写入数据
	for (unsigned int i = 0; i < length; i++)
	{
		MyPoint mp(i + 1, (type == ROW) ? *POINTELEM(data, cols, index, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index));
		ret = WriteFile(handle_write, &mp, sizeof(MyPoint), &len, NULL); //子进程读了后，父进程才可以继续写入管道
		if (!ret)
		{
			_ERROR("绘制图像", "写入数据点失败");
			return;
		}
	}
	//关闭管道
	if (!CloseHandle(handle_write)) 
	{
		_ERROR("绘制图像", "关闭管道失败");
		return;
	}
}


void MatBase::Plot(unsigned int type, unsigned int index1, unsigned int index2, const char* tittle)
{
	//判断异常
	if (!OnInit)
	{
		_ERROR("绘制图像", "未初始化");
		return;
	}
	if (!PathFileExistsA(PATH_PLOT))
	{
		char szBuffer[64];
		sprintf(szBuffer, "找不到画图模块 %s", PATH_PLOT);
		_ERROR("绘制图像", szBuffer);
		return;
	}
	if ((type == ROW && GetMax(index1, index2) > rows) || (type == COL && GetMax(index1, index2) > cols))
	{
		_ERROR("绘制图像", "坐标超限");
		return;
	}
	//采样压缩
	unsigned int length0 = (type == ROW) ? cols : rows;
	unsigned int sumb = GetBitNum(length0);
	unsigned int ite = 1;
	if (sumb > 3)
		ite = pow(10, sumb - 3.0);
	unsigned int length = length0 / ite;
	//进程通信管道
	SECURITY_ATTRIBUTES saAttr;
	saAttr.nLength = sizeof(SECURITY_ATTRIBUTES);
	saAttr.bInheritHandle = TRUE;
	saAttr.lpSecurityDescriptor = NULL;
	HANDLE handle_read;
	HANDLE handle_write;
	bool ret = CreatePipe(&handle_read, &handle_write, &saAttr, 0);
	if (!ret)
	{
		_ERROR("绘制图像", "创建进程通信管道失败");
		return;
	}
	//启动绘图模块进程
	char cmdline[] = PATH_PLOT;
	PROCESS_INFORMATION piProcInfo;
	ZeroMemory(&piProcInfo, sizeof(PROCESS_INFORMATION));
	STARTUPINFO si;
	ZeroMemory(&si, sizeof(STARTUPINFO));
	si.cb = sizeof(STARTUPINFO);
	si.hStdInput = handle_read; //把管道的读句柄传给子进程
	si.dwFlags |= STARTF_USESTDHANDLES;
	ret = CreateProcess(NULL, cmdline, NULL, NULL, TRUE, CREATE_NEW_CONSOLE, NULL, NULL, &si, &piProcInfo);
	if (!ret)
	{
		_ERROR("绘制图像", "启动绘图模块进程失败");
		return;
	}
	else
	{
		CloseHandle(piProcInfo.hProcess);
		CloseHandle(piProcInfo.hThread);
		CloseHandle(handle_read);
	}
	//写入文件头
	FigHead fh;
	memcpy(fh.szTitle, tittle, (strlen(tittle) + 1) * sizeof(char));
	fh.length = length;
	DWORD len;
	ret = WriteFile(handle_write, &fh, sizeof(FigHead), &len, NULL); //子进程读了后，父进程才可以继续写入管道
	if (!ret)
	{
		_ERROR("绘制图像", "写入文件头失败");
		return;
	}
	//写入数据
	for (unsigned int i = 0; i < length; i++)
	{
		MyPoint mp((type == ROW) ? *POINTELEM(data, cols, index1, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index1), (type == ROW) ? *POINTELEM(data, cols, index2, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index2));
		ret = WriteFile(handle_write, &mp, sizeof(MyPoint), &len, NULL); //子进程读了后，父进程才可以继续写入管道
		if (!ret)
		{
			_ERROR("绘制图像", "写入数据点失败");
			return;
		}
	}
	//关闭管道
	if (!CloseHandle(handle_write))
	{
		_ERROR("绘制图像", "关闭管道失败");
		return;
	}
}

void MatBase::PlotD(unsigned int type, unsigned int index, unsigned int index1, unsigned int index2, const char* tittle)
{
	//判断异常
	if (!OnInit)
	{
		_ERROR("绘制图像", "未初始化");
		return;
	}
	if (!PathFileExistsA(PATH_PLOT))
	{
		char szBuffer[64];
		sprintf(szBuffer, "找不到画图模块 %s", PATH_PLOT);
		_ERROR("绘制图像", szBuffer);
		return;
	}
	if ((type == ROW && GetMax3(index, index1, index2) > rows) || (type == COL && GetMax3(index, index1, index2) > cols))
	{
		_ERROR("绘制图像", "坐标超限");
		return;
	}
	//采样压缩
	unsigned int length0 = (type == ROW) ? cols : rows;
	unsigned int sumb = GetBitNum(length0);
	unsigned int ite = 1;
	if (sumb > 3)
		ite = pow(10, sumb - 3.0);
	unsigned int length = length0 / ite;
	//进程通信管道
	SECURITY_ATTRIBUTES saAttr;
	saAttr.nLength = sizeof(SECURITY_ATTRIBUTES);
	saAttr.bInheritHandle = TRUE;
	saAttr.lpSecurityDescriptor = NULL;
	HANDLE handle_read;
	HANDLE handle_write;
	bool ret = CreatePipe(&handle_read, &handle_write, &saAttr, 0);
	if (!ret)
	{
		_ERROR("绘制图像", "创建进程通信管道失败");
		return;
	}
	//启动绘图模块进程
	char cmdline[] = PATH_PLOT;
	PROCESS_INFORMATION piProcInfo;
	ZeroMemory(&piProcInfo, sizeof(PROCESS_INFORMATION));
	STARTUPINFO si;
	ZeroMemory(&si, sizeof(STARTUPINFO));
	si.cb = sizeof(STARTUPINFO);
	si.hStdInput = handle_read; //把管道的读句柄传给子进程
	si.dwFlags |= STARTF_USESTDHANDLES;
	ret = CreateProcess(NULL, cmdline, NULL, NULL, TRUE, CREATE_NEW_CONSOLE, NULL, NULL, &si, &piProcInfo);
	if (!ret)
	{
		_ERROR("绘制图像", "启动绘图模块进程失败");
		return;
	}
	else
	{
		CloseHandle(piProcInfo.hProcess);
		CloseHandle(piProcInfo.hThread);
		CloseHandle(handle_read);
	}
	//写入文件头
	FigHead fh;
	memcpy(fh.szTitle, tittle, (strlen(tittle) + 1) * sizeof(char));
	fh.length = 2 * length;
	DWORD len;
	ret = WriteFile(handle_write, &fh, sizeof(FigHead), &len, NULL); //子进程读了后，父进程才可以继续写入管道
	if (!ret)
	{
		_ERROR("绘制图像", "写入文件头失败");
		return;
	}
	//写入数据
	for (unsigned int i = 0; i < length; i++)
	{
		MyPoint mp((type == ROW) ? *POINTELEM(data, cols, index, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index), (type == ROW) ? *POINTELEM(data, cols, index1, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index1));
		ret = WriteFile(handle_write, &mp, sizeof(MyPoint), &len, NULL); //子进程读了后，父进程才可以继续写入管道
		if (!ret)
		{
			_ERROR("绘制图像", "写入数据点失败");
			return;
		}
	}
	for (unsigned int i = 0; i < length; i++)
	{
		MyPoint mp((type == ROW) ? *POINTELEM(data, cols, index, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index), (type == ROW) ? *POINTELEM(data, cols, index2, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index2));
		ret = WriteFile(handle_write, &mp, sizeof(MyPoint), &len, NULL); //子进程读了后，父进程才可以继续写入管道
		if (!ret)
		{
			_ERROR("绘制图像", "写入数据点失败");
			return;
		}
	}
	//关闭管道
	if (!CloseHandle(handle_write))
	{
		_ERROR("绘制图像", "关闭管道失败");
		return;
	}
}

void MatBase::PlotD(unsigned int type, unsigned int index, unsigned int index1, unsigned int index2, unsigned int index3, const char* tittle)
{
	//判断异常
	if (!OnInit)
	{
		_ERROR("绘制图像", "未初始化");
		return;
	}
	if (!PathFileExistsA(PATH_PLOT))
	{
		char szBuffer[64];
		sprintf(szBuffer, "找不到画图模块 %s", PATH_PLOT);
		_ERROR("绘制图像", szBuffer);
		return;
	}
	if ((type == ROW && GetMax3(index, index1, index2) > rows) || (type == COL && GetMax3(index, index1, index2) > cols))
	{
		_ERROR("绘制图像", "坐标超限");
		return;
	}
	//采样压缩
	unsigned int length0 = (type == ROW) ? cols : rows;
	unsigned int sumb = GetBitNum(length0);
	unsigned int ite = 1;
	if (sumb > 3)
		ite = pow(10, sumb - 3.0);
	unsigned int length = length0 / ite;
	//进程通信管道
	SECURITY_ATTRIBUTES saAttr;
	saAttr.nLength = sizeof(SECURITY_ATTRIBUTES);
	saAttr.bInheritHandle = TRUE;
	saAttr.lpSecurityDescriptor = NULL;
	HANDLE handle_read;
	HANDLE handle_write;
	bool ret = CreatePipe(&handle_read, &handle_write, &saAttr, 0);
	if (!ret)
	{
		_ERROR("绘制图像", "创建进程通信管道失败");
		return;
	}
	//启动绘图模块进程
	char cmdline[] = PATH_PLOT;
	PROCESS_INFORMATION piProcInfo;
	ZeroMemory(&piProcInfo, sizeof(PROCESS_INFORMATION));
	STARTUPINFO si;
	ZeroMemory(&si, sizeof(STARTUPINFO));
	si.cb = sizeof(STARTUPINFO);
	si.hStdInput = handle_read; //把管道的读句柄传给子进程
	si.dwFlags |= STARTF_USESTDHANDLES;
	ret = CreateProcess(NULL, cmdline, NULL, NULL, TRUE, CREATE_NEW_CONSOLE, NULL, NULL, &si, &piProcInfo);
	if (!ret)
	{
		_ERROR("绘制图像", "启动绘图模块进程失败");
		return;
	}
	else
	{
		CloseHandle(piProcInfo.hProcess);
		CloseHandle(piProcInfo.hThread);
		CloseHandle(handle_read);
	}
	//写入文件头
	FigHead fh;
	memcpy(fh.szTitle, tittle, (strlen(tittle) + 1) * sizeof(char));
	fh.length = 3 * length;
	DWORD len;
	ret = WriteFile(handle_write, &fh, sizeof(FigHead), &len, NULL); //子进程读了后，父进程才可以继续写入管道
	if (!ret)
	{
		_ERROR("绘制图像", "写入文件头失败");
		return;
	}
	//写入数据
	for (unsigned int i = 0; i < length; i++)
	{
		MyPoint mp((type == ROW) ? *POINTELEM(data, cols, index, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index), (type == ROW) ? *POINTELEM(data, cols, index1, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index1));
		ret = WriteFile(handle_write, &mp, sizeof(MyPoint), &len, NULL); //子进程读了后，父进程才可以继续写入管道
		if (!ret)
		{
			_ERROR("绘制图像", "写入数据点失败");
			return;
		}
	}
	for (unsigned int i = 0; i < length; i++)
	{
		MyPoint mp((type == ROW) ? *POINTELEM(data, cols, index, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index), (type == ROW) ? *POINTELEM(data, cols, index2, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index2));
		ret = WriteFile(handle_write, &mp, sizeof(MyPoint), &len, NULL); //子进程读了后，父进程才可以继续写入管道
		if (!ret)
		{
			_ERROR("绘制图像", "写入数据点失败");
			return;
		}
	}
	for (unsigned int i = 0; i < length; i++)
	{
		MyPoint mp((type == ROW) ? *POINTELEM(data, cols, index, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index), (type == ROW) ? *POINTELEM(data, cols, index3, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index3));
		ret = WriteFile(handle_write, &mp, sizeof(MyPoint), &len, NULL); //子进程读了后，父进程才可以继续写入管道
		if (!ret)
		{
			_ERROR("绘制图像", "写入数据点失败");
			return;
		}
	}
	//关闭管道
	if (!CloseHandle(handle_write))
	{
		_ERROR("绘制图像", "关闭管道失败");
		return;
	}
}

void MatBase::Plot3(unsigned int type, unsigned int index1, unsigned int index2, unsigned int index3, const char* tittle)
{
	//判断异常
	if (!OnInit)
	{
		_ERROR("绘制图像", "未初始化");
		return;
	}
	if (!PathFileExistsA(PATH_PLOT3))
	{
		char szBuffer[64];
		sprintf(szBuffer, "找不到画图模块 %s", PATH_PLOT3);
		_ERROR("绘制图像", szBuffer);
		return;
	}
	if ((type == ROW && GetMax3(index1, index2, index3) > rows) || (type == COL && GetMax3(index1, index2, index3) > cols))
	{
		_ERROR("绘制图像", "坐标超限");
		return;
	}
	//采样压缩
	unsigned int length0 = (type == ROW) ? cols : rows;
	unsigned int sumb = GetBitNum(length0);
	unsigned int ite = 1;
	if (sumb > 2)
		ite = pow(10, sumb - 2.0);
	unsigned int length = length0 / ite;
	//进程通信管道
	SECURITY_ATTRIBUTES saAttr;
	saAttr.nLength = sizeof(SECURITY_ATTRIBUTES);
	saAttr.bInheritHandle = TRUE;
	saAttr.lpSecurityDescriptor = NULL;
	HANDLE handle_read;
	HANDLE handle_write;
	bool ret = CreatePipe(&handle_read, &handle_write, &saAttr, 0);
	if (!ret)
	{
		_ERROR("绘制图像", "创建进程通信管道失败");
		return;
	}
	//启动绘图模块进程
	char cmdline[] = PATH_PLOT3;
	PROCESS_INFORMATION piProcInfo;
	ZeroMemory(&piProcInfo, sizeof(PROCESS_INFORMATION));
	STARTUPINFO si;
	ZeroMemory(&si, sizeof(STARTUPINFO));
	si.cb = sizeof(STARTUPINFO);
	si.hStdInput = handle_read; //把管道的读句柄传给子进程
	si.dwFlags |= STARTF_USESTDHANDLES;
	ret = CreateProcess(NULL, cmdline, NULL, NULL, TRUE, CREATE_NEW_CONSOLE, NULL, NULL, &si, &piProcInfo);
	if (!ret)
	{
		_ERROR("绘制图像", "启动绘图模块进程失败");
		return;
	}
	else
	{
		CloseHandle(piProcInfo.hProcess);
		CloseHandle(piProcInfo.hThread);
		CloseHandle(handle_read);
	}
	//写入文件头
	FigHead fh;
	memcpy(fh.szTitle, tittle, (strlen(tittle) + 1) * sizeof(char));
	fh.length = length;
	DWORD len;
	ret = WriteFile(handle_write, &fh, sizeof(FigHead), &len, NULL); //子进程读了后，父进程才可以继续写入管道
	if (!ret)
	{
		_ERROR("绘制图像", "写入文件头失败");
		return;
	}
	//写入数据
	for (unsigned int i = 0; i < length; i++)
	{
		MyPoint3 mp((type == ROW) ? *POINTELEM(data, cols, index1, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index1), (type == ROW) ? *POINTELEM(data, cols, index2, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index2), (type == ROW) ? *POINTELEM(data, cols, index3, 1 + i * ite) : *POINTELEM(data, cols, 1 + i * ite, index3));
		ret = WriteFile(handle_write, &mp, sizeof(MyPoint3), &len, NULL); //子进程读了后，父进程才可以继续写入管道
		if (!ret)
		{
			_ERROR("绘制图像", "写入数据点失败");
			return;
		}
	}
	//关闭管道
	if (!CloseHandle(handle_write))
	{
		_ERROR("绘制图像", "关闭管道失败");
		return;
	}
}

MatBase operator+(const MatBase& m1, MatBase m2) 
{
	InitTemp(m2);
	if (m1.rows != m2.rows || m1.cols != m2.cols) 
	{
		_ERROR("矩阵相加","行列数不匹配");
		return MatBase();
	}
	MatBase ans;
	ans = m1;
	for (int i = 0;i < ans.rows * ans.cols;i++)
		ans.data[i] += m2.data[i];
	DelTemp(m2);
	ans.OnInit = 2;
	return ans;
}
MatBase operator-(const MatBase& m1, MatBase m2)
{
	InitTemp(m2);
	if (m1.rows != m2.rows || m1.cols != m2.cols)
	{
		_ERROR("矩阵相减", "行列数不匹配");
		return MatBase();
	}
	MatBase ans;
	ans = m1;
	for (int i = 0;i < ans.rows * ans.cols;i++)
		ans.data[i] -= m2.data[i];
	DelTemp(m2);
	ans.OnInit = 2;
	return ans;
}
MatBase operator+(const real& val, const MatBase& mat) 
{
	MatBase ans;
	ans = mat;
	for (int i = 0;i < ans.rows * ans.cols;i++)
		ans.data[i] += val;
	ans.OnInit = 2;
	return ans;
}
MatBase operator+(const MatBase& mat, const real& val) 
{
	MatBase ans;
	ans = mat;
	for (int i = 0;i < ans.rows * ans.cols;i++)
		ans.data[i] += val;
	ans.OnInit = 2;
	return ans;
}
MatBase operator-(const real& val, const MatBase& mat) 
{
	MatBase ans;
	ans = mat;
	for (int i = 0;i < ans.rows * ans.cols;i++)
		ans.data[i] = val - ans.data[i];
	ans.OnInit = 2;
	return ans;
}
MatBase operator-(const MatBase& mat, const real& val) 
{
	MatBase ans;
	ans = mat;
	for (int i = 0;i < ans.rows * ans.cols;i++)
		ans.data[i] -= val;
	ans.OnInit = 2;
	return ans;
}
MatBase operator*(const real& val, const MatBase& mat)
{
	MatBase ans;
	ans = mat;
	for (int i = 0;i < ans.rows * ans.cols;i++)
		ans.data[i] *= val;
	ans.OnInit = 2;
	return ans;
}
MatBase operator*(const MatBase& mat, const real& val) 
{
	MatBase ans;
	ans = mat;
	for (int i = 0;i < ans.rows * ans.cols;i++)
		ans.data[i] *= val;
	ans.OnInit = 2;
	return ans;
}
MatBase operator/(const MatBase& mat, const real& val) 
{
	MatBase ans;
	ans = mat;
	for (int i = 0;i < ans.rows * ans.cols;i++)
		ans.data[i] /= val;
	ans.OnInit = 2;
	return ans;
}
MatBase& MatBase::operator=(MatBase mat)
{
	if (OnInit && data)
		free(data);
	OnInit = 1;
	rows = mat.rows;
	cols = mat.cols;
	if (mat.OnInit == 1) //复制永久对象
	{ 
		data = (real*)malloc(mat.rows * mat.cols * sizeof(real));
		memcpy(data, mat.data, rows * cols * sizeof(real));
	}
	else if (mat.OnInit == 2)  //传递地址，销毁临时对象
	{
		data = mat.data;
		//related = true;
	}
	else
		_ERROR("对象赋值", "赋值右值未初始化");
	return *this;
}
bool operator==(const MatBase& m1, const MatBase& m2) 
{
	if (m1.rows != m2.rows || m1.cols != m2.cols)
	{
		_ERROR("矩阵判断大小", "行列数不匹配");
		return false;
	}
	for (unsigned int i = 0;i < m1.rows * m1.cols;i++)
		if (!Equal(m1.data[i], m2.data[i]))
			return false;
	return true;
}
MatBase operator*(MatBase m1, MatBase m2) 
{
	InitTemp(m1);
	InitTemp(m2);
	if (m1.cols != m2.rows) 
	{
		_ERROR("矩阵相乘","行列数不匹配");
		return MatBase();
	}
	MatBase ans(m1.rows, m2.cols);
	ans.SetAll(0);
	for (int x = 1;x <= ans.rows;x++)
		for (int y = 1;y <= ans.cols;y++)
			for (int n = 1;n <= m1.cols;n++)
				*POINTELEM(ans.data, ans.cols, x, y) += (*POINTELEM(m1.data, m1.cols, x, n)) * (*POINTELEM(m2.data, m2.cols, n, y));
	DelTemp(m1);
	DelTemp(m2);
	ans.OnInit = 2;
	return ans;
}
MatBase MatBase::MultipAll(const MatBase& m1, MatBase m2)
{
	InitTemp(m2);
	if (m1.rows != m2.rows || m1.cols != m2.cols) 
	{
		_ERROR("矩阵对应元素依次相乘","矩阵行列数不匹配");
		return MatBase();
	}
	MatBase ans;
	ans = m1;
	for (int i = 0;i < ans.rows * ans.cols;i++)
		ans.data[i] *= m2.data[i];
	DelTemp(m2);
	ans.OnInit = 2;
	return ans;
}
/*实现Vector中方法*/
Vector::Vector() {}
Vector::Vector(unsigned int _rows, unsigned int _cols, int _OnInit) :MatBase(_rows, _cols, _OnInit)
{
	if (_rows > 1 && _cols > 1)
	{
		_ERROR("创建向量", "行数大于1或列数大于1");
		OnInit = 0;
		rows = 0;
		cols = 0;
		free(data);
		data = NULL;
		//related = false;
	}
}
Vector::Vector(real low, real dis, real up) :MatBase(1, round((up - low) / dis) + 1, true)
{
	if ((up - low) * dis <= 0)
	{
		_ERROR("初始化均匀采样向量", "上下界限大小和间距不匹配");
		return;
	}
	unsigned int index = 1;
	if (dis > 0)
		for (real val = low;val <= up + EqualLimit;val += dis)
		{
			*POINTELEMV(data, index) = val;
			index++;
		}
	else
		for (real val = low;val >= up - EqualLimit;val += dis)
		{
			*POINTELEMV(data, index) = val;
			index++;
		}
}
Vector::Vector(const MatBase& vec)
{
	OnInit = 1;
	rows = vec.rows;
	cols = vec.cols;
	if (vec.OnInit == 1) //复制永久对象
	{
		data = (real*)malloc(rows * cols * sizeof(real));
		memcpy(data, vec.data, rows * cols * sizeof(real));
	}
	else if (vec.OnInit == 2)  //传递地址，销毁临时对象
	{
		data = vec.data;
	}
	else
		_ERROR("对象赋值", "赋值右值未初始化");
}
int Vector::Type() const
{
	if (rows == 1)
		return 1;
	if (cols == 1)
		return 2;
	return 0;
}
unsigned int Vector::Length() const
{

	unsigned int n;
	if (Type() == 1)
		n = cols;
	else
		n = rows;
	return n;

}
real Vector::Norm() const
{
	real sum = 0;
	for (int i = 0;i < Length();i++)
		sum += data[i] * data[i];
	return sqrt(sum);
}
Vector Vector::GetTrans() const
{
	Vector v(cols, rows, 0);
	v.data = data;
	v.OnInit = OnInit;
	return v;
}
bool Vector::SetTrans()
{
	int t = rows;
	rows = cols;
	cols = t;
	return true;
}
void Vector::SetUni(real low, real dis, real up) 
{
	if ((up - low) * dis <= 0) 
	{
		_ERROR("获得均匀采样向量","上下界限大小和间距不匹配");
		return;
	}
	if (OnInit) 
		free(data);
	OnInit = 1;
	rows = 1;
	cols = round((up - low) / dis) + 1;
	data = (real*)malloc(rows * cols * sizeof(real));
	unsigned int index = 1;
	if (dis > 0) 
		for (real val = low;val <= up + EqualLimit;val += dis)
		{
			*POINTELEMV(data, index) = val;
			index++;
		}
	else 
		for (real val = low;val >= up - EqualLimit;val += dis)
		{
			*POINTELEMV(data, index) = val;
			index++;
		}
}
Vector Vector::GetPart(unsigned int st, unsigned int ed) const
{
	if (st > ed) 
	{
		_ERROR("截取向量", "开始坐标大于终点坐标");
		return Vector();
	}
	unsigned int _row, _cols;
	if (Type() == 1) 
	{
		_row = 1;
		_cols = ed - st + 1;
	}
	else 
	{
		_row = ed - st + 1;
		_cols = 1;
	}
	Vector ans(_row, _cols);
	memcpy(ans.data, data + st - 1, (ed - st + 1) * sizeof(real));
	ans.OnInit = 2;
	return ans;
}
Vector Vector::GetPart(Vector index) const
{
	InitTemp(index);
	if (!index.OnInit)
	{
		_ERROR("截取向量", "下标向量未初始化");
		return Vector();
	}
	unsigned int _row, _cols;
	if (Type() == 1)
	{
		_row = 1;
		_cols = index.Length();
	}
	else
	{
		_row = index.Length();
		_cols = 1;
	}
	Vector ans(_row, _cols);
	for (int i = 1;i <= ans.Length();i++)
		*POINTELEMV(ans.data, i) = *POINTELEMV(data, (int)*POINTELEMV(index.data, i));
	DelTemp(index);
	ans.OnInit = 2;
	return ans;
}
Vector& Vector::SetRotateX(real rad)
{
	if (Length() != 3) 
	{
		_ERROR("向量旋转", "向量维数不为3");
		return *this;
	}
	real crad = cos(rad);
	real srad = sin(rad);
	Matrix m(3,3);
	m[1][1] = 1; m[1][2] = 0; m[1][3] = 0;
	m[2][1] = 0; m[2][2] = crad; m[2][3] = srad;
	m[3][1] = 0; m[3][2] = -srad; m[3][3] = crad;
	*this = m * *this;
	m.Delete();
	return *this;
}
Vector& Vector::SetRotateY(real rad)
{
	if (Length() != 3)
	{
		_ERROR("向量旋转", "向量维数不为3");
		return *this;
	}
	real crad = cos(rad);
	real srad = sin(rad);
	Matrix m(3, 3);
	m[1][1] = crad; m[1][2] = 0; m[1][3] = -srad;
	m[2][1] = 0; m[2][2] = 1; m[2][3] = 0;
	m[3][1] = srad; m[3][2] = 0; m[3][3] = crad;
	*this = m * *this;
	m.Delete();
	return *this;
}
Vector& Vector::SetRotateZ(real rad)
{
	if (Length() != 3)
	{
		_ERROR("向量旋转", "向量维数不为3");
		return *this;
	}
	real crad = cos(rad);
	real srad = sin(rad);
	Matrix m(3, 3);
	m[1][1] = crad; m[1][2] = srad; m[1][3] = 0;
	m[2][1] = -srad; m[2][2] = crad; m[2][3] = 0;
	m[3][1] = 0; m[3][2] = 0; m[3][3] = 1;
	*this = m * *this;
	m.Delete();
	return *this;
}
Vector Vector::GetRotateX(real rad)
{
	if (Length() != 3)
	{
		_ERROR("向量旋转", "向量维数不为3");
		return Vector();
	}
	real crad = cos(rad);
	real srad = sin(rad);
	Matrix m(3, 3);
	m[1][1] = 1; m[1][2] = 0; m[1][3] = 0;
	m[2][1] = 0; m[2][2] = crad; m[2][3] = srad;
	m[3][1] = 0; m[3][2] = -srad; m[3][3] = crad;
	Vector ans;
	ans = *this;
	ans = m * ans;
	m.Delete();
	ans.OnInit = 2;
	return ans;
}
Vector Vector::GetRotateY(real rad)
{
	if (Length() != 3)
	{
		_ERROR("向量旋转", "向量维数不为3");
		return Vector();
	}
	real crad = cos(rad);
	real srad = sin(rad);
	Matrix m(3, 3);
	m[1][1] = crad; m[1][2] = 0; m[1][3] = -srad;
	m[2][1] = 0; m[2][2] = 1; m[2][3] = 0;
	m[3][1] = srad; m[3][2] = 0; m[3][3] = crad;
	Vector ans;
	ans = *this;
	ans = m * ans;
	m.Delete();
	ans.OnInit = 2;
	return ans;
}
Vector Vector::GetRotateZ(real rad)
{
	if (Length() != 3)
	{
		_ERROR("向量旋转", "向量维数不为3");
		return Vector();
	}
	real crad = cos(rad);
	real srad = sin(rad);
	Matrix m(3, 3);
	m[1][1] = crad; m[1][2] = srad; m[1][3] = 0;
	m[2][1] = -srad; m[2][2] = crad; m[2][3] = 0;
	m[3][1] = 0; m[3][2] = 0; m[3][3] = 1;
	Vector ans;
	ans = *this;
	ans = m * ans;
	m.Delete();
	ans.OnInit = 2;
	return ans;
}

void Vector::GetArray(real *array, unsigned int size)
{
	memcpy(array, data, size * sizeof(real));
}

void Vector::FromArray(real* array, unsigned int size, int type)
{
	if (OnInit)
		free(data);
	OnInit = 1;
	if (type == 1)
	{
		rows = 1;
		cols = size;
	}
	else
	{
		rows = size;
		cols = 1;
	}
	data = (real*)malloc(size * sizeof(real));
	memcpy(data, array, size * sizeof(real));
}

Vector& Vector::operator=(MatBase vec)
{
	if (OnInit && data)
		free(data);
	OnInit = 1;
	rows = vec.rows;
	cols = vec.cols;
	if (vec.OnInit == 1) //复制永久对象
	{
		data = (real*)malloc(rows * cols * sizeof(real));
		memcpy(data, vec.data, rows * cols * sizeof(real));
	}
	else if (vec.OnInit == 2)  //传递地址，销毁临时对象
	{
		data = vec.data;
		//related = true;
	}
	else
		_ERROR("对象赋值", "赋值右值未初始化");
	return *this;
}

real Vector::Dot(Vector v1, Vector v2)
{
	InitTemp(v1);
	InitTemp(v2);
	if (v1.Length() != v2.Length())
	{
		_ERROR("向量点乘", "向量长度不匹配");
		return 0;
	}
	real sum = 0;
	for (int i = 0;i < v1.Length();i++)
		sum += v1.data[i] * v2.data[i];
	DelTemp(v1);
	DelTemp(v2);
	return sum;
}
Vector Vector::Cross(Vector v1, Vector v2, int type) 
{
	InitTemp(v1);
	InitTemp(v2);
	if (v1.Length() != 3 || v2.Length() != 3) 
	{
		_ERROR("向量叉乘", "向量维度不为3");
		return Vector();
	}
	int rows, cols;
	if (type == 1) 
	{
		rows = 1;
		cols = v1.Length();
	}
	else
	{
		rows = v1.Length();
		cols = 1;
	}
	Vector ans(rows, cols);
	ans[1] = v1[2] * v2[3] - v1[3] * v2[2];
	ans[2] = v1[3] * v2[1] - v1[1] * v2[3];
	ans[3] = v1[1] * v2[2] - v1[2] * v2[1];
	DelTemp(v1);
	DelTemp(v2);
	ans.OnInit = 2;
	return ans;
}
Vector Vector::Cat(Vector v1, Vector v2 , int type)
{
	InitTemp(v1);
	InitTemp(v2);
	unsigned int rows, cols;
	if (!type) 
		if (v1.Type() != v2.Type()) 
		{
			_ERROR("向量拼接", "向量类型不匹配");
			return Vector();
		}
		else 
		{
			type = v1.Type();
		}
	if (type == 1) //行向量
	{
		rows = 1;
		cols = v1.Length() + v2.Length();
	}
	else //列向量
	{
		rows = v1.Length() + v2.Length();
		cols = 1;
	}
	Vector ans(rows, cols);
	memcpy(ans.data, v1.data, v1.Length() * sizeof(real));
	memcpy(POINTELEMV(ans.data, v1.Length() + 1), v2.data, v2.Length() * sizeof(real));
	DelTemp(v1);
	DelTemp(v2);
	ans.OnInit = 2;
	return ans;
}
real& Vector::operator[](const unsigned int index)
{
	return data[index - 1];
}

/*实现Matrix中方法*/
Matrix::Matrix() 
{
	OnLU = 0;
	L = NULL;
	U = NULL;
}
Matrix::Matrix(unsigned int _rows, unsigned int _cols, int _OnInit) :MatBase(_rows, _cols, _OnInit)
{
	OnLU = 0;
	L = NULL;
	U = NULL;
}
Matrix::Matrix(const MatBase& mat)
{
	OnInit = 1;
	rows = mat.rows;
	cols = mat.cols;
	if (mat.OnInit == 1) //复制永久对象
	{
		data = (real*)malloc(mat.rows * mat.cols * sizeof(real));
		memcpy(data, mat.data, rows * cols * sizeof(real));
	}
	else if (mat.OnInit == 2)  //传递地址，销毁临时对象
	{
		data = mat.data;
	}
	else
		_ERROR("对象赋值", "赋值右值未初始化");
}
void Matrix::Init(unsigned int _rows, unsigned int _cols, int _OnInit)
{
	OnInit = _OnInit;
	rows = _rows;
	cols = _cols;
	OnLU = 0;
	L = NULL;
	U = NULL;
	if (_OnInit)
		data = (real*)malloc(rows * cols * sizeof(real));
	//related = false;
}
void Matrix::Delete() 
{
	if (OnLU)
	{
		OnLU = 0;
		if (L)
			free(L);
		if (U)
			free(U);
	}
	if (OnInit)
	{
		OnInit = 0;
		rows = 0;
		cols = 0;
		if (data)
			free(data);
		data = NULL;
	}
}
void Matrix::DeleteLU() 
{
	if (OnLU)
	{
		OnLU = 0;
		if (L)
			free(L);
		if (U)
			free(U);
	}
}
bool Matrix::InitLU()
{
	if (rows != cols)
	{
		_ERROR("初始化LU矩阵", "矩阵不是方阵");
		return false;
	}
	if (OnLU)
	{
		OnLU = 0;
		free(L);
		free(U);
	}
	OnLU = rows;
	L = (real*)malloc(OnLU * OnLU * sizeof(real));
	U = (real*)malloc(OnLU * OnLU * sizeof(real));
	//memset(L, 0, OnLU * OnLU * sizeof(real));
	//memset(U, 0, OnLU * OnLU * sizeof(real));
	for (int j = 1;j <= OnLU;j++)
	{
		*POINTELEM(U, OnLU, 1, j) = *POINTELEM(data, OnLU, 1, j);
		*POINTELEM(L, OnLU, j, j) = 1;
	}
	if (fabs(*POINTELEM(U, OnLU, 1, 1)) < NanLimit)
		*POINTELEM(U, OnLU, 1, 1) = NanLimit;
	for (int i = 2;i <= OnLU;i++)
		*POINTELEM(L, OnLU, i, 1) = *POINTELEM(data, OnLU, i, 1) / *POINTELEM(U, OnLU, 1, 1);
	for (int k = 2;k <= OnLU;k++) 
	{
		for (int j = k;j <= OnLU;j++)
		{
			real sum = 0;
			for (int m = 1;m <= k - 1;m++)
				sum += *POINTELEM(L, OnLU, k, m) * *POINTELEM(U, OnLU, m, j);
			*POINTELEM(U, OnLU, k, j) = *POINTELEM(data, OnLU, k, j) - sum;
		}
		for (int i = k + 1;i <= OnLU;i++)
		{
			real sum = 0;
			for (int m = 1;m <= k - 1;m++)
				sum += *POINTELEM(L, OnLU, i, m) * *POINTELEM(U, OnLU, m, k);
			if (fabs(*POINTELEM(U, OnLU, k, k)) < NanLimit)
				*POINTELEM(U, OnLU, k, k) = NanLimit;
			*POINTELEM(L, OnLU, i, k) = (*POINTELEM(data, OnLU, i, k) - sum) / *POINTELEM(U, OnLU, k, k);
		}
	}
	return true;
}
real* Matrix::PointRow(unsigned int x)
{
	return POINTROW(data, cols, x);
}
Vector Matrix::GetRow(unsigned int x, bool copy) const
{
	if (x > rows) 
	{
		_ERROR("矩阵取行","坐标超限");
		return Vector();
	}
	if (copy) 
	{
		Vector v(1, cols);
		memcpy(v.data, POINTROW(data, cols, x), cols * sizeof(real));
		v.OnInit = 2;
		return v;
	}
	else 
	{
		Vector v(1, cols, 0);
		v.data = POINTROW(data, cols, x);
		return v;
	}
}

Vector Matrix::GetCol(unsigned int y) const
{
	if (y > cols)
	{
		_ERROR("矩阵取列", "坐标超限");
		return Vector();
	}
	Vector v(rows, 1);
	for (int i = 1;i <= rows;i++)
		v[i] = *POINTELEM(data, cols, i, y);
	v.OnInit = 2;
	return v;
}
bool Matrix::SetRow(unsigned int x, Vector val)
{
	InitTemp(val);
	if (cols != val.Length())
	{
		_ERROR("矩阵按行赋值", "矩阵列数与向量长度不匹配");
		return false;
	}
	if (x > rows)
	{
		_ERROR("矩阵按行赋值", "坐标超限");
		return false;
	}
	memcpy(POINTROW(data, cols, x), val.data, cols * sizeof(real));
	DelTemp(val);
	return true;
}
bool Matrix::SetCol(unsigned int y, Vector val) 
{
	InitTemp(val);
	if (rows != val.Length())
	{
		_ERROR("矩阵按列赋值", "矩阵行数与向量长度不匹配");
		return false;
	}
	if (y > cols)
	{
		_ERROR("矩阵按列赋值", "坐标超限");
		return false;
	}
	for (int i = 1;i <= rows;i++)
		*POINTELEM(data, cols, i, y) = *POINTELEMV(val.data, i);
	DelTemp(val);
	return true;
}
bool Matrix::SetEye()
{
	if (rows != cols)
	{
		_ERROR("初始化LU矩阵", "矩阵不是方阵");
		return false;
	}
	memset(data, 0, rows * cols * sizeof(real));
	for (int i = 1;i <= rows;i++)
		*POINTELEM(data, cols, i, i) = 1;
	return true;
}
Vector Matrix::LinSolve(Vector b) const
{
	InitTemp(b);
	if (!OnLU)
	{
		_ERROR("解线性方程组", "LU矩阵未初始化");
		return Vector();
	}
	if (OnLU != b.Length())
	{
		_ERROR("解线性方程组", "向量维度不匹配");
		return Vector();
	}
	if (b.Type() == 2) 
	{
		Vector y(OnLU, 1);
		*POINTELEMV(y.data, 1) = *POINTELEMV(b.data, 1);
		for (int k = 2;k <= OnLU;k++)
		{
			real sum = 0;
			for (int m = 1;m <= k - 1;m++)
				sum += *POINTELEM(L, OnLU, k, m) * *POINTELEMV(y.data, m);
			*POINTELEMV(y.data, k) = *POINTELEMV(b.data, k) - sum;
		}

		Vector x(OnLU, 1);
		*POINTELEMV(x.data, OnLU) = *POINTELEMV(y.data, OnLU) / *POINTELEM(U, OnLU, OnLU, OnLU);
		for (int k = OnLU - 1;k >= 1;k--)
		{
			real sum = 0;
			for (int m = k + 1;m <= OnLU;m++)
				sum += *POINTELEM(U, OnLU, k, m) * *POINTELEMV(x.data, m);
			*POINTELEMV(x.data, k) = (*POINTELEMV(y.data, k) - sum) / *POINTELEM(U, OnLU, k, k);
		}
		y.Delete();
		DelTemp(b);
		x.OnInit = 2;
		return x;
	}
	else 
	{
		Vector y(1, OnLU);
		*POINTELEMV(y.data, 1) = *POINTELEMV(b.data, 1) / *POINTELEM(U, OnLU, 1, 1);
		for (int k = 2;k <= OnLU;k++)
		{
			real sum = 0;
			for (int m = 1;m <= k - 1;m++)
				sum += *POINTELEM(U, OnLU, m, k) * *POINTELEMV(y.data, m);
			*POINTELEMV(y.data, k) = (*POINTELEMV(b.data, k) - sum) / *POINTELEM(U, OnLU, k, k);
		}

		Vector x(1, OnLU);
		*POINTELEMV(x.data, OnLU) = *POINTELEMV(y.data, OnLU);
		for (int k = OnLU - 1;k >= 1;k--)
		{
			real sum = 0;
			for (int m = k + 1;m <= OnLU;m++)
				sum += *POINTELEM(L, OnLU, m, k) * *POINTELEMV(x.data, m);
			*POINTELEMV(x.data, k) = *POINTELEMV(y.data, k) - sum;
		}
		y.Delete();
		DelTemp(b);
		x.OnInit = 2;
		return x;
	}
}
Matrix Matrix::GetInv(bool OnInit)
{
	if (rows != cols)
	{
		_ERROR("获得逆矩阵", "矩阵不是方阵");
		return Matrix();
	}
	if (OnInit)
		InitLU();
	Matrix eye(rows, cols);
	eye.SetEye();
	Matrix ans(rows, cols);
	ans = DiviRight(eye, *this);
	eye.Delete();
	ans.OnInit = 2;
	return ans;
}

bool Matrix::SetInv(bool OnInit)
{
	if (rows != cols)
	{
		_ERROR("矩阵求逆", "矩阵不是方阵");
		return false;
	}
	if (OnInit)
		InitLU();
	Matrix eye(rows, cols);
	eye.SetEye();
	Matrix ans(rows, cols);
	ans = DiviRight(eye, *this);
	memcpy(data, ans.data, rows * cols * sizeof(real));
	ans.Delete();
	eye.Delete();
	return true;
}

Matrix Matrix::GetPart(int type, unsigned int st, unsigned int ed) const
{
	if (st > ed)
	{
		_ERROR("截取矩阵", "开始坐标大于终点坐标");
		return Matrix();
	}
	unsigned int _row, _cols;
	if (type == ROW) //按行截取
	{
		_row = ed - st + 1;
		_cols = cols;
	}
	else //按列截取
	{
		_row = rows;
		_cols = ed - st + 1;
	}
	Matrix ans(_row, _cols);
	if (type == ROW) //按行截取
		for (int i = 1;i <= ed - st + 1;i++)
			ans.SetRow(i, this->GetRow(st + i - 1));
	else
		for (int i = 1;i <= ed - st + 1;i++)
			ans.SetCol(i, this->GetCol(st + i - 1));
	ans.OnInit = 2;
	return ans;
}
/*
Vector Vector::GetPart(unsigned int st, unsigned int ed) const
{
	if (st > ed)
	{
		_ERROR("截取向量", "开始坐标大于终点坐标");
		return Vector();
	}
	unsigned int _row, _cols;
	if (Type() == 1)
	{
		_row = 1;
		_cols = ed - st + 1;
	}
	else
	{
		_row = ed - st + 1;
		_cols = 1;
	}
	Vector ans(_row, _cols);
	memcpy(ans.data, data + st - 1, (ed - st + 1) * sizeof(real));
	ans.OnInit = 2;
	return ans;
}
*/

Vector Matrix::operator[](const unsigned int index)
{
	Vector v(1, cols , 0);
	v.data = POINTROW(data, cols, index);
	return v;
}
Matrix& Matrix::operator=(MatBase mat)
{
	if (OnInit && data)
		free(data);
	OnInit = 1;
	rows = mat.rows;
	cols = mat.cols;
	if (mat.OnInit == 1) //复制永久对象
	{
		data = (real*)malloc(mat.rows * mat.cols * sizeof(real));
		memcpy(data, mat.data, rows * cols * sizeof(real));
	}
	else if (mat.OnInit == 2)  //传递地址，销毁临时对象
	{
		data = mat.data;
		//related = true;
	}
	else
		_ERROR("对象赋值", "赋值右值未初始化");
	return *this;
}

Matrix Matrix::DiviLeft(Matrix A, Matrix B)
{
	InitTemp(A);
	InitTemp(B);
	if (A.cols != B.rows) 
	{
		_ERROR("矩阵左除","维度不匹配");
		return Matrix();
	}
	Matrix X(B.rows, B.cols);
	for (int y = 1;y <= X.cols;y++) 
		X.SetCol(y, A.LinSolve(B.GetCol(y)));
	DelTemp(A);
	DelTemp(B);
	X.OnInit = 2;
	return X;
}
Matrix Matrix::DiviRight(Matrix B, Matrix A)
{
	InitTemp(B);
	InitTemp(A);
	if (A.rows != B.cols)
	{
		_ERROR("矩阵右除", "维度不匹配");
		return Matrix();
	}
	Matrix X(B.rows, B.cols);
	for (int x = 1;x <= X.rows;x++)
		X.SetRow(x, A.LinSolve(B.GetRow(x)));
	DelTemp(B);
	DelTemp(A);
	X.OnInit = 2;
	return X;
}
Matrix Matrix::GetRMX(real rad)
{
	real crad = cos(rad);
	real srad = sin(rad);
	Matrix m(3, 3);
	m[1][1] = 1; m[1][2] = 0; m[1][3] = 0;
	m[2][1] = 0; m[2][2] = crad; m[2][3] = srad;
	m[3][1] = 0; m[3][2] = -srad; m[3][3] = crad;
	m.OnInit = 2;
	return m;
}
Matrix Matrix::GetRMY(real rad)
{
	real crad = cos(rad);
	real srad = sin(rad);
	Matrix m(3, 3);
	m[1][1] = crad; m[1][2] = 0; m[1][3] = -srad;
	m[2][1] = 0; m[2][2] = 1; m[2][3] = 0;
	m[3][1] = srad; m[3][2] = 0; m[3][3] = crad;
	m.OnInit = 2;
	return m;
}
Matrix Matrix::GetRMZ(real rad)
{
	real crad = cos(rad);
	real srad = sin(rad);
	Matrix m(3, 3);
	m[1][1] = crad; m[1][2] = srad; m[1][3] = 0;
	m[2][1] = -srad; m[2][2] = crad; m[2][3] = 0;
	m[3][1] = 0; m[3][2] = 0; m[3][3] = 1;
	m.OnInit = 2;
	return m;
}




/*实现Legendre中方法*/
Legendre::Legendre(unsigned int _N) :BASE(_N, _N)
{
	OnInit = false;
	N = _N;
}
void Legendre::Delete()
{
	BASE.Delete();
}
void Legendre::InitBASE()
{
	BASE.SetZeros();
	BASE[1][1] = 1; //n=1 1
	BASE[2][2] = 1; //n=2 s
	Vector p1(1, N), p2(1, N), p(1, N);
	for (int i = 3;i <= N;i++)
	{
		p2 = (i - 2.) / (i - 1.) * BASE.GetRow(i - 2);
		p1 = (2. * i - 3.) / (i - 1.) * BASE.GetRow(i - 1);
		p.SetZeros();
		memcpy(p.data + 1, p1.data, (i - 1) * sizeof(real));
		p = p - p2;
		BASE.SetRow(i, p);
	}
	OnInit = true;
}
Vector Legendre::AddCoe(const Vector& C)
{
	if (!OnInit)
	{
		_ERROR("勒让德正交基系数叠加", "正交基未初始化");
		return Vector();
	}
	if (C.Length() != N)
	{
		_ERROR("勒让德正交基系数叠加", "系数向量维度不匹配");
		return Vector();
	}
	Vector sum(1, N);
	sum.SetZeros();
	for (int i = 1;i <= N;i++)
		sum = *POINTELEMV(C.data, i) * BASE.GetRow(i) + sum;
	return sum;
}
Vector Legendre::Fitting(const Vector& s, const Vector& f)
{
	Vector C(1, N), phi(1, N), mul(1, N);
	for (int i = 0;i < N;i++)
	{
		phi = PolyValue(BASE[i + 1], s);
		mul = MatBase::MultipAll(phi, f);
		C[i + 1] = Integral(s, mul) * (2. * i + 1.) / 2.;
	}
	return C;
}
/*函数定义*/
real PolyValue(Vector A, real x)
{
	InitTemp(A);
	real val = 0;
	for (int i = A.Length();i >= 1;i--) 
	{
		val *= x;
		val += *POINTELEMV(A.data, i);
	}
	DelTemp(A);
	return val;
}
Vector PolyValue(Vector A, Vector x)
{
	InitTemp(A);
	InitTemp(x);
	Vector val(x.rows, x.cols);
	val.SetZeros();
	for (int i = A.Length();i >= 1;i--)
	{
		val = MatBase::MultipAll(val, x);
		val = val + *POINTELEMV(A.data, i);
	}
	DelTemp(A);
	DelTemp(x);
	val.OnInit = 2;
	return val;
}
real Integral(Vector x, Vector f)
{
	InitTemp(x);
	InitTemp(f);
	if (x.Length() != f.Length()) 
	{
		_ERROR("数值积分","变量维度不匹配");
		return 0;
	}
	real val = 0;
	for (int i = 1;i < x.Length();i++)
		val += (*POINTELEMV(f.data, i + 1) + *POINTELEMV(f.data, i)) * (*POINTELEMV(x.data, i + 1) - *POINTELEMV(x.data, i)) / 2;
	DelTemp(x);
	DelTemp(f);
	return val;
}

unsigned int GetBitNum(double x)
{
	if (x < 0)
		x *= -1;
	unsigned int sumb = 0;
	if (x > 1)
		while (x > 10)
		{
			x /= 10;
			sumb++;
		}
	else
		while (x < 1)
		{
			x *= 10;
			sumb--;
		}
	return sumb;
}